Ejercicio Final - Series Temporales
Javier Herráez Albarrán
22/4/2022
Carga de librerías y datos
library(ggplot2)
library(plotly)
library(tseries)
library(MASS)
library(forecast)
library(lmtest)
library(caschrono)
library(lubridate)
library(timeDate)
library(tsoutliers)
library(dplyr)
# Se limpia el espacio de trabajo:
rm(list = ls())
# Se asigna el directorio donde están los datos
# setwd("C:/Users/Javier/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")
setwd("C:/Users/jherraez/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")
datos <- read.csv("_11_CACERES.csv")Dividimos en train y en set y se convierten los datos en un objeto de tipo time series para poder aplicar sobre ellos las funciones de R.
datos$FECHA = as.Date(datos$FECHA,format='%d/%m/%Y')
datos.train <- subset(datos, FECHA<=as.Date('01/12/2018',format='%d/%m/%Y'))
datos.test <- subset(datos, FECHA>as.Date('01/12/2018',format='%d/%m/%Y'))
datos.train.ts <- as.ts(datos.train$TARGET, frequency = 12)
datos.test.ts <- as.ts(datos.test$TARGET, frequency=12)Representación gráfica de los datos
names(datos.train)## [1] "FECHA" "TARGET"graficoInicial <- ggplot(aes(x= FECHA, y = TARGET), data = datos.train) +
geom_line(color = '#d84519', size = 1) +
xlab('FECHA') + ylab('Matriculaciones')
ggplotly(graficoInicial)Estacionalidad
Estacionalidad en varianza
Evaluamos si es necesario transformar la serie para hacerla estacionaria en varianza.
box_cox <- boxcox(TARGET ~ FECHA,
data = datos.train,
lambda = c(0, 0.5, 1))
lambda <- box_cox$x[which.max(box_cox$y)]
lambda## [1] 0.1616162Como lambda toma un valor cercano a cero, tomamos logaritmos.
datos.train$log_target=log(datos.train$TARGET)
datos.test$log_target=log(datos.test$TARGET)
datos.train.ts <- as.ts(datos.train$log_target, frequency = 12)
datos.test.ts <- as.ts(datos.test$log_target, frequency=12)Estacionalidad en media
Probamos a hacer el test de Dickey-Fuller aunque no confiamos demasiado en su resultado, ya que es poco exigente.
adf.test(datos.train.ts)##
## Augmented Dickey-Fuller Test
##
## data: datos.train.ts
## Dickey-Fuller = -1.6659, Lag order = 5, p-value = 0.7159
## alternative hypothesis: stationaryTodo parece indicar que va a ser necesaria una diferencia regular, pero nosotros ajustamos directamente sin nada.
Ajuste de un modelo SARIMA a la serie
Observamos los gráficos f.a.s. y la f.a.p. para realizar el primer ajuste.
acf(datos.train.ts, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(datos.train.ts, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
Lo habitual es empezar con una estructura AR(1). La fuerte correlación de orden 1 en la f.a.p parece reforzar esta decisión, por lo tanto, se propone un SARIMA(1,0,0)x(0,0,0)12 + MU
ajuste1 <- Arima(datos.train.ts,
order = c(1,0,0),
seasonal = list(order = c(0,0,0), period = 12),
method = "ML")
ajuste1## Series: datos.train.ts
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.9059 5.5001
## s.e. 0.0303 0.2199
##
## sigma^2 = 0.09146: log likelihood = -42.67
## AIC=91.34 AICc=91.47 BIC=101.11coeftest(ajuste1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.905855 0.030348 29.849 < 2.2e-16 ***
## intercept 5.500117 0.219941 25.007 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajuste1) # debajo de 0.8 todo, no correlaciones## ar1 intercept
## ar1 1.00000000 -0.02691084
## intercept -0.02691084 1.00000000Box.test.2(residuals(ajuste1),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box") # queremos valores >0.05 en todos## Retard p-value
## [1,] 6 6.49e-06
## [2,] 12 5.80e-07
## [3,] 18 2.54e-06
## [4,] 24 1.10e-07
## [5,] 30 6.80e-07
## [6,] 36 8.00e-08
## [7,] 42 4.30e-07
## [8,] 48 3.00e-07acf(ajuste1$residuals, lag.max = 48, xlab = "Retardo", # aquí vemos que sólo se salen en multiplos de 12 y que hay patrón de simetría # para la practica es medio verdad, decidimos hacer lo otro
main= "Función de autocorrelación simple")
pacf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
# Se propone un SARIMA(1,0,1)x(0,0,0)12 + MUajuste2 <- Arima(datos.train.ts,
order = c(1,0,1),
seasonal = list(order = c(0,0,0), period = 12),
method = "ML")
coeftest(ajuste2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.9947456 0.0054474 182.610 < 2.2e-16 ***
## ma1 -0.7254036 0.0434252 -16.705 < 2.2e-16 ***
## intercept 5.5238319 0.5351063 10.323 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajuste2)## ar1 ma1 intercept
## ar1 1.0000000 -0.27585575 -0.01079400
## ma1 -0.2758557 1.00000000 -0.01006758
## intercept -0.0107940 -0.01006758 1.00000000Box.test.2(residuals(ajuste2),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box") # test de jung-box## Retard p-value
## [1,] 6 0.69073487
## [2,] 12 0.22973902
## [3,] 18 0.21893065
## [4,] 24 0.03960481
## [5,] 30 0.08546994
## [6,] 36 0.03580800
## [7,] 42 0.04708855
## [8,] 48 0.06789183acf(ajuste2$residuals, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(ajuste2$residuals, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
# Se propone un SARIMA(1,0,1)x(1,0,0)12 + MUajuste3 <- Arima(datos.train.ts,
order = c(1,0,1),
seasonal = list(order = c(1,0,0), period = 12),
method = "ML")
coeftest(ajuste3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.9922742 0.0072296 137.2515 < 2.2e-16 ***
## ma1 -0.7393339 0.0458966 -16.1087 < 2.2e-16 ***
## sar1 0.2426026 0.0734512 3.3029 0.0009569 ***
## intercept 5.5194141 0.4959124 11.1298 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Se propone un SARIMA(0,1,1)x(1,0,0)12 + MUajuste4 <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
method = "ML")
coeftest(ajuste4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.743005 0.044523 -16.6881 < 2.2e-16 ***
## sar1 0.238164 0.072912 3.2664 0.001089 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajuste4)## ma1 sar1
## ma1 1.0000000 -0.1216682
## sar1 -0.1216682 1.0000000Box.test.2(residuals(ajuste4),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box") # test de jung-box## Retard p-value
## [1,] 6 0.7123801
## [2,] 12 0.9055981
## [3,] 18 0.8138609
## [4,] 24 0.5663195
## [5,] 30 0.7499105
## [6,] 36 0.6405039
## [7,] 42 0.6366659
## [8,] 48 0.7534805calculoExplicativasCalendarioCaceres <- function(variableFecha, domingoYFestivosJuntos){
#######################################
# Creación de todas las fechas #
#######################################
if (month(max(variableFecha)) %in% c(1,3,5,7,8,10,12)) {
diasHastaFinMes <- 30
} else if (month(max(variableFecha)) %in% c(4,6,9,11)) {
diasHastaFinMes <- 29
} else if (year(max(variableFecha))%%4==0) {
diasHastaFinMes <- 28
} else {diasHastaFinMes <- 27}
todasLasFechas <- data.frame(fechas=seq(min(variableFecha),
max(variableFecha)+diasHastaFinMes,
by="days"))
#######################################
# Cálculo de la Semana Santa #
#######################################
domingoResurrecion <- as.Date(Easter(year(min(variableFecha)):year(max(variableFecha))))
lunesPascua <- domingoResurrecion + 1
sabadoSanto <- domingoResurrecion - 1
viernesSanto <- domingoResurrecion - 2
juevesSanto <- domingoResurrecion - 3
# Se unen y ordenan todos los días que forman la Semana Santa
semanaSanta <- sort(c(juevesSanto, viernesSanto, sabadoSanto, domingoResurrecion, lunesPascua))
# Se pone en formato data.frame y se añade un indicador
semanaSanta <- data.frame(fechas=semanaSanta, semanaSanta=rep(1,length(semanaSanta)))
# Se añaden a la tabla maestra de fechas
todasLasFechas_2 <- merge(x = todasLasFechas, y = semanaSanta, by = "fechas", all.x = TRUE)
# Se reemplazan los NAs por 0
todasLasFechas_2$semanaSanta[is.na(todasLasFechas_2$semanaSanta)] <- 0
######################################
# Cálculo de la variable dt #
######################################
# 1. Definición de festivos:
############################
calendario <- todasLasFechas
calendario$diaSemana <- as.factor(wday(calendario$fecha))
calendario$diaMes <- as.factor(day(calendario$fecha))
calendario$mes <- as.factor(month(calendario$fecha))
calendario$anyo <- as.factor(year(calendario$fecha))
calendario$p_01ene <- ifelse(calendario$diaMes==1 & calendario$mes==1, 1, 0)
calendario$p_06ene <- ifelse(calendario$diaMes==6 & calendario$mes==1, 1, 0)
calendario$p_19mar <- ifelse(calendario$diaMes==19 & calendario$mes==3, 1, 0)
calendario$p_01may <- ifelse(calendario$diaMes==1 & calendario$mes==5, 1, 0)
calendario$p_15ago <- ifelse(calendario$diaMes==15 & calendario$mes==8, 1, 0)
calendario$p_12oct <- ifelse(calendario$diaMes==12 & calendario$mes==10,1, 0)
calendario$p_01nov <- ifelse(calendario$diaMes==1 & calendario$mes==11, 1 ,0)
calendario$p_06dic <- ifelse(calendario$diaMes==6 & calendario$mes==12, 1 ,0)
calendario$p_08dic <- ifelse(calendario$diaMes==8 & calendario$mes==12, 1 ,0)
calendario$p_25dic <- ifelse(calendario$diaMes==25 & calendario$mes==12, 1 ,0)
# Fiestas regionales: día de Extremadura
calendario$p_08sep <- ifelse(calendario$diaMes==8 & calendario$mes==9, 1 ,0)
# Fiestas locales: San Jorge
calendario$p_23abr <- ifelse(calendario$diaMes==23 & calendario$mes==4, 1 ,0)
calendario$festivo <- rowSums(subset(calendario, select=p_01ene:p_25dic))
# La definición de la variable dt varía según la opción domingoYFestivosJuntos.
if (domingoYFestivosJuntos==0){
calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
# Días laborables: todos menos sábados y domingos
calendario$laborable <- 1-calendario$sabado-calendario$domingo
} else {
calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
# Domingo=1 si domingo=1 o festivo=1
calendario$domingo <- ifelse(calendario$domingo==1 | calendario$festivo==1, 1 ,0)
# Días laborables: todos menos sábados y domingos/festivos
calendario$laborable <- 1-calendario$sabado-calendario$domingo
}
# 2. Definición de variable dt:
###############################
# Se filtran las columnas de inter?s y se añade la Semana Santa
calendario_2 <- calendario[, c("fechas", "mes", "anyo", "sabado", "domingo", "laborable", "festivo")]
todasLasFechasFinal <- merge(x = todasLasFechas_2, y = calendario_2,
by = "fechas", all.x = TRUE)
# Agregamos la serie a nivel a?o-mes
calendarioAnyoMes <- aggregate(todasLasFechasFinal[,c("sabado","domingo",
"laborable", "semanaSanta", "festivo")],
by=list(mes=todasLasFechasFinal$mes,
anyo=todasLasFechasFinal$anyo),
"sum")
# Se calcula la variable dt:
calendarioAnyoMes$dt <- calendarioAnyoMes$laborable-(5/2)*(calendarioAnyoMes$sabado+calendarioAnyoMes$domingo)
######################################
# Cálculo de años bisiestos #
######################################
calendarioAnyoMes$anyoNum <- as.numeric(levels(calendarioAnyoMes$anyo))[calendarioAnyoMes$anyo]
calendarioAnyoMes$bisiesto <- ifelse(calendarioAnyoMes$mes==2 &(calendarioAnyoMes$anyoNum %% 4)==0, 1 ,0)
#######################################################
# Tabla final con explicativas de calendario #
#######################################################
if (domingoYFestivosJuntos==0){
explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto", "festivo")])
} else {
explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto")])
}
return(explicativasCalendario)
}explicativasCalendarioTrain <- calculoExplicativasCalendarioCaceres(datos.train$FECHA, domingoYFestivosJuntos=0)
calendarioTrain <-
as.matrix(
explicativasCalendarioTrain[,c("semanaSanta", "dt", "bisiesto", "festivo")]
)ajuste4conCalendario <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioTrain,
method = "ML")
coeftest(ajuste4conCalendario)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7199518 0.0483371 -14.8944 < 2.2e-16 ***
## sar1 0.1968911 0.0733748 2.6834 0.007289 **
## semanaSanta -0.0115610 0.0118108 -0.9788 0.327654
## dt 0.0150418 0.0055963 2.6878 0.007192 **
## bisiesto -0.0227074 0.1049431 -0.2164 0.828693
## festivo -0.0614426 0.0211852 -2.9003 0.003729 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# decidimos quitar año bisiestocalendarioTrain <-
as.matrix(
explicativasCalendarioTrain[,c("semanaSanta", "dt", "festivo")]
)
# Se podría incluir tambi?n la variable festivo
ajuste5conCalendario <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioTrain,
method = "ML")
coeftest(ajuste5conCalendario)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7202627 0.0482284 -14.9344 < 2.2e-16 ***
## sar1 0.1975910 0.0732627 2.6970 0.006996 **
## semanaSanta -0.0113950 0.0117906 -0.9664 0.333821
## dt 0.0150369 0.0055984 2.6859 0.007233 **
## festivo -0.0607887 0.0209877 -2.8964 0.003775 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajuste5conCalendario)## ma1 sar1 semanaSanta dt festivo
## ma1 1.000000000 -0.1395694285 -0.02512203 -0.0060496156 -0.07048090
## sar1 -0.139569429 1.0000000000 -0.07245826 -0.0004302792 -0.03481278
## semanaSanta -0.025122031 -0.0724582600 1.00000000 0.0491736995 0.16027535
## dt -0.006049616 -0.0004302792 0.04917370 1.0000000000 0.00428523
## festivo -0.070480901 -0.0348127821 0.16027535 0.0042852305 1.00000000Box.test.2(residuals(ajuste5conCalendario),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")## Retard p-value
## [1,] 6 0.7984990
## [2,] 12 0.9116931
## [3,] 18 0.8911888
## [4,] 24 0.7504270
## [5,] 30 0.9289244
## [6,] 36 0.8699718
## [7,] 42 0.8579565
## [8,] 48 0.9146346# [7]. Identificaci?n de outliers
# Esto identifica outliers
listaOutliersTrain <- locate.outliers(ajuste5conCalendario$residuals,
pars = coefs2poly(ajuste5conCalendario),
types = c("AO", "LS", "TC"),
cval=3)
listaOutliersTrain$abststat=abs(listaOutliersTrain$tstat)
# Cruzamos con la tabla original para obtener la fecha
datos.train$ind <- as.numeric(rownames(datos.train))
listaOutliersTrainFecha <- merge(listaOutliersTrain, datos.train[,c("ind", "FECHA")], by = "ind")
arrange(listaOutliersTrainFecha,desc(listaOutliersTrainFecha$abststat))## ind type coefhat tstat abststat FECHA
## 1 99 LS -0.4456849 -3.135621 3.135621 2011-03-01
## 2 126 TC -0.5369467 -3.119897 3.119897 2013-06-01
## 3 88 TC 0.5255352 3.053591 3.053591 2010-04-01outliersTrain <- outliers(c( "TC","LS","AO"), c(88, 99, 126))
outliersVariablesTrain <- outliers.effects(outliersTrain, length(ajuste5conCalendario$residuals))
calendarioMasOutliersTrain <- as.matrix(cbind(calendarioTrain, outliersVariablesTrain))ajuste5conCalendarioYOutliers <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioMasOutliersTrain,
method = "ML")
coeftest(ajuste5conCalendarioYOutliers)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7705343 0.0443086 -17.3902 < 2.2e-16 ***
## sar1 0.2043685 0.0742688 2.7517 0.005928 **
## semanaSanta -0.0133834 0.0115593 -1.1578 0.246947
## dt 0.0138329 0.0055408 2.4966 0.012541 *
## festivo -0.0611110 0.0203547 -3.0023 0.002679 **
## TC88 0.4155507 0.1833772 2.2661 0.023445 *
## LS99 -0.3817126 0.1443567 -2.6442 0.008188 **
## AO126 -0.4052362 0.2049407 -1.9773 0.048004 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.arma(ajuste5conCalendarioYOutliers)## ma1 sar1 semanaSanta dt festivo
## ma1 1.00000000 -0.099074246 -0.03840997 0.016422714 -0.033633401
## sar1 -0.09907425 1.000000000 -0.05525868 -0.011715920 -0.030682112
## semanaSanta -0.03840997 -0.055258676 1.00000000 0.054744004 0.152545681
## dt 0.01642271 -0.011715920 0.05474400 1.000000000 0.009458039
## festivo -0.03363340 -0.030682112 0.15254568 0.009458039 1.000000000
## TC88 0.12700224 -0.003696872 -0.12348004 -0.011857780 0.053345083
## LS99 0.10033436 0.160271312 -0.07843167 -0.026932065 0.038140304
## AO126 0.16167463 0.024750256 0.01876841 0.137038092 0.043949322
## TC88 LS99 AO126
## ma1 0.127002243 0.10033436 0.16167463
## sar1 -0.003696872 0.16027131 0.02475026
## semanaSanta -0.123480045 -0.07843167 0.01876841
## dt -0.011857780 -0.02693207 0.13703809
## festivo 0.053345083 0.03814030 0.04394932
## TC88 1.000000000 0.25658620 0.02030734
## LS99 0.256586201 1.00000000 0.02203153
## AO126 0.020307344 0.02203153 1.00000000Box.test.2(residuals(ajuste5conCalendarioYOutliers),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")## Retard p-value
## [1,] 6 0.9386016
## [2,] 12 0.9353550
## [3,] 18 0.8785053
## [4,] 24 0.8525545
## [5,] 30 0.9582087
## [6,] 36 0.9668843
## [7,] 42 0.9358235
## [8,] 48 0.9096749[7]. Predicción
explicativasCalendarioTest <- calculoExplicativasCalendarioCaceres(datos.test$FECHA,domingoYFestivosJuntos=0)
calendarioTest <- as.matrix(explicativasCalendarioTest[,c("semanaSanta", "dt", "festivo")])# Nos valen los ?ltimos 12 datos de la variable (12 1's)
# para el valor futuro de la variable
outliersVariablesTest <- outliersVariablesTrain[181:192,]
calendarioMasOutliersTest <- as.matrix(cbind(calendarioTest, outliersVariablesTest))
calendarioMasOutliersTest## semanaSanta dt festivo TC88 LS99 AO126
## [1,] 0 3.0 2 3.927514e-15 1 0
## [2,] 0 0.0 0 2.749260e-15 1 0
## [3,] 0 -4.0 1 1.924482e-15 1 0
## [4,] 5 2.0 0 1.347137e-15 1 0
## [5,] 0 3.0 1 9.429961e-16 1 0
## [6,] 0 -5.0 0 6.600972e-16 1 0
## [7,] 0 3.0 0 4.620681e-16 1 0
## [8,] 0 -0.5 1 3.234477e-16 1 0
## [9,] 0 -1.5 0 2.264134e-16 1 0
## [10,] 0 3.0 1 1.584893e-16 1 0
## [11,] 0 -1.5 1 1.109425e-16 1 0
## [12,] 0 -0.5 3 7.765978e-17 1 0prediccion <- as.data.frame(predict(ajuste5conCalendarioYOutliers,
n.ahead=12,
newxreg=calendarioMasOutliersTest))
#L?mites de confianza al 95%
U <- exp(prediccion$pred + 2*prediccion$se)
L <- exp(prediccion$pred - 2*prediccion$se)
datos.pred <- data.frame(FECHA = datos.test$FECHA,
Prediccion = exp(prediccion$pred + 0.5 * prediccion$se^2),
LimSup = U,
LimInf = L)datos.real.pred <- merge(datos, datos.pred, by = "FECHA", all.x = T)
datos.real.pred$REAL <- datos.real.pred$TARGET
datos.real.pred$TARGET[datos.real.pred$FECHA>as.Date('01/12/2018',format='%d/%m/%Y')] <- NA
grafico <- ggplot(data = datos.real.pred) +
geom_line(aes(x= FECHA, y = REAL), color = 'blue',
alpha = 0.4, size = 0.8) +
geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
alpha = 0.8, size = 0.8) +
geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
size = 1) +
geom_line(aes(x= FECHA, y = LimSup), color = 'orange',
size = 1) +
geom_line(aes(x= FECHA, y = LimInf), color = 'orange',
size = 1) +
xlab('FECHA') + ylab('Matriculaciones')
ggplotly(grafico)